from IPython.display import Image
from IPython.core.display import HTML
Image(url= "https://img.techpowerup.org/201119/logo147.png", width=900, height=300)
Abbreviations: \ GMP : Green Mountain Power ( distributor) \ PRO : Prosumer
The decision variables' set then is defined as : $X_t \doteq \{ x^{RL}_{t}, x^{RB}_{t}, x^{RE}_{t}, x^{RG}_{t}, x^{BG}_{t}, x^{EG}_{t}, x^{BL}_{t}, x^{EL}_{t}, x^{GL}_{t}, x^{GB}_{t}, x^{GE}_{t}, x^{peak}_{k}\}$
s.t. \begin{align} x^{GL}_t + \eta^d (x^{BL}_t + x^{EL}_t) + x^{RL}_t \geq& L_t & \forall t \label{c_load_2} \tag{2}\\ x^{BL}_t + x^{BG}_t \leq& b_t^B & \forall t \label{c_bat_state_2} \tag{3}\\ x^{EL}_t + x^{EG}_t \leq& b_t^E A_t & \forall t \label{c_ev_state_2} \tag{4}\\ x^{GE}_t + x^{RE}_t \leq& E^{max} A_t & \forall t \label{c_ev_2_state_2} \tag{5}\\ x^{BL}_t + x^{EL}_t + x^{BG}_t + x^{EG}_t \leq& U^d & \forall t \label{c_dch_max_2} \tag{6}\\ x^{GB}_t + x^{GE}_t + x^{RB}_t + x^{RE}_t \leq& U^c & \forall t \label{c_ch_max_2} \tag{7}\\ B^{min} \leq b_t^B \leq& B^{max} & \forall t \label{c_b_bat_max_2} \tag{8}\\ E^{min} \leq b_t^E \leq& E^{max} & \forall t \label{c_b_ev_max_2} \tag{9}\\ x^{RL}_t + x^{RB}_t + x^{RE}_t + x^{RG}_t \leq& R_t & \forall t \label{c_gen_cap_2} \tag{10}\\ b_t^B + \eta^c (x^{GB}_t + x^{RB}_t) - ( x^{BL}_t + x^{BG}_t) =& b_{t+1}^B & \forall t \label{c_bat_next_state_2} \tag{11}\\ b_t^E + \eta^c (x^{GE}_t + x^{RE}_t) - ( x^{EL}_t + x^{EG}_t) - \frac{V_t}{\eta^d} =& b_{t+1}^E & \forall t \label{c_ev_next_state_2} \tag{12}\\ x^{GL}_{t} + x^{GB}_{t} + x^{GE}_{t} \leq& x_{k}^{peak} & \forall t \in k, k\label{c_max_2} \tag{13}\\ x^{RL}_{t}, x^{RB}_{t}, x^{RE}_{t}, x^{RG}_{t}, x^{BG}_{t}, x^{EG}_{t}, x^{BL}_{t}, x^{EL}_{t}, x^{GL}_{t}, x^{GB}_{t}, x^{GE}_{t}, x^{peak}_{k}, b_t^B, b_t^E \geq& 0 & \forall t \label{c_domain_2} \tag{14} \end{align}
Objective function \ref{obj_2} minimizes electricity consumption cost. Constraint \ref{c_load_2} makes sure load is covered. Constraints \ref{c_bat_state_2} to \ref{c_ev_2_state_2} state that electricity taken out of battery or EV must be firstly less than the available energy stored in them and secondly, electricity can be taken from or sent to EV, when it is available. Constraints \ref{c_dch_max_2} and \ref{c_ch_max_2} indicate the charging and discharging limits. Constraints \ref{c_b_bat_max_2} and \ref{c_b_ev_max_2} make sure the state of battery and EV are in the allowed range of stored energy in them. Constraint \ref{c_gen_cap_2} checks electricity sent from solar panel would be less than or equal to the generated electricity and constraints \ref{c_bat_next_state_2} and \ref{c_ev_next_state_2} calculate the state of battery and EV in the next episode. Constraint \ref{c_max_2} calculates the maximum electricity taken from grid in month $k$. Finally, constraint \ref{c_domain_2} specifies the domain of the decision variables.
s.t. constraints 2 t0 12 and \begin{align} x^{RG}_t + x^{BG}_t + x^{EG}_t \leq& (U^d + R_t) I^{trade} & \forall t\label{c_max_3} \tag{16}\\ x^{RL}_{t}, x^{RB}_{t}, x^{RE}_{t}, x^{RG}_{t}, x^{BG}_{t}, x^{EG}_{t}, x^{BL}_{t}, x^{EL}_{t}, x^{GL}_{t}, x^{GB}_{t}, x^{GE}_{t}, b_t^B, b_t^E \geq& 0 & \forall t \label{c_domain_3} \tag{17} \end{align}
Objective function \ref{obj_3} minimizes electricity consumption cost. Constraint \ref{c_max_3} makes sure the prosumer can sell electricity to grid if trade is allowed and, constraint \ref{c_domain_3} specifies the domain of the decision variables.
from model import *
import pandas as pd
import numpy as np
import plotly.express as px
import plotly.graph_objects as go
from plotly.subplots import make_subplots
scenarios = ["s0","s1","s2","s3","s4","s5"]
pd.options.display.max_columns = 500
pd.options.display.max_rows = 200
pd.options.display.float_format = '{:,.0f}'.format
import gurobipy as gp
from gurobipy import GRB
import pickle
import numpy_financial as npf
pro_results, pro_dres, pro_key_results, pro_npvs, pro_npvs2 = pro_ctrl(run=False,house_num=4)
gmp_results, gmp_dres, gmp_key_results, gmp_npvs, gmp_npvs2= gmp_ctrl(run=False,house_num=4)
fig = px.line(pro_dres[(pro_dres.Rate=='Time-of-Use') & (pro_dres.Sell==True) & (pro_dres.Scenario=='s1')], x='Time (Month-Day/Hour:Minute)', y=["P (GMP)","P (Flat)","P (ToU)", "P (ToU for EV)"], animation_frame= 'Scenario', title="Electricity price", labels={"value":"$", "P (GMP)": "P (RTLMP)"}, range_y=[0,0.3])
newnames = {'P (GMP)':'RTLMP', 'P (Flat)': 'Flat', 'P (ToU)': 'ToU', 'P (ToU for EV)': 'ToU for EV'}
fig.for_each_trace(lambda t: t.update(name = newnames[t.name],
legendgroup = newnames[t.name],
hovertemplate = t.hovertemplate.replace(t.name, newnames[t.name])
)
)
fig.show()
gmp_dres= gmp_dres.replace(['s0','s1','s2','s3','s4','s5','Time-of-Use', "Time-of-Use EV & Flat", "flat"],['s0: Status Quo', 's1: Battery', 's2: Battery+PV', 's3: EV', 's4: EV+PV', 's5: EV+PV+Battery', "TOU", "TOU EV & Flat", "Flat"])
tmp = gmp_dres[gmp_dres.Scenario == 's3: EV']
# Create figure with secondary y-axis
fig = make_subplots(specs=[[{"secondary_y": True}]])
# Add traces
fig.add_trace(
go.Scatter(x=tmp['Time (Month-Day/Hour:Minute)'], y=tmp["Flow from Grid"], name="Flow from Grid"),
secondary_y=False,
)
fig.add_trace(
go.Scatter(x=tmp['Time (Month-Day/Hour:Minute)'], y=tmp["P_GMP"], name="RTLMP"),
secondary_y=True,
)
# fig.update_layout(
# title_text="Double Y Axis Example"
# )
fig.update_xaxes(title_text="Time (Month-Day/Hour:Minute)")
# Set y-axes titles
fig.update_yaxes(title_text="<b>kWh</b> (flow)", secondary_y=False)
fig.update_yaxes(title_text="<b>$/kWh</b> (rate)", secondary_y=True)
fig.show()
print('The dirstributor controls')
for i in range(0,6):
db = gmp_results[gmp_results.Scenario==f"s{i}"]
fig = go.Figure(data=[go.Sankey(
node = dict(
pad = 15,
thickness = 20,
line = dict(color = "black", width = 0.5),
label = ["Grid (G)", "Solar PV (R)", "ESS (B)", "Load (L)", "EV (E)"]
),
link = dict(
source = [0, 0, 1, 1, 1, 2, 2, 0, 1, 4, 4],
target = [2, 3, 0, 2, 3, 0, 3, 4, 4, 0, 3],
value = [db.loc[:, "x^GB"].sum(), db.loc[:, "x^GL"].sum(), db.loc[:, "x^RG"].sum(), db.loc[:, "x^RB"].sum(), db.loc[:, "x^RL"].sum(), db.loc[:, "x^BG"].sum(), db.loc[:, "x^BL"].sum(), db.loc[:, "x^GE"].sum(), db.loc[:, "x^RE"].sum(), db.loc[:, "x^EG"].sum(), db.loc[:, "x^EL"].sum()]
))])
fig.update_layout(title_text=f"Scenario {i}", font_size=10)
fig.show()
pro_dres = pro_dres.replace(['s0','s1','s2','s3','s4','s5','Time-of-Use', "Time-of-Use EV & Flat", "flat"],['s0: Status Quo', 's1: Battery', 's2: Battery+PV', 's3: EV', 's4: EV+PV', 's5: EV+PV+Battery', "TOU", "TOU EV & Flat", "Flat"])
fig = px.line(pro_dres[(pro_dres.Rate=='TOU') & (pro_dres.Sell==True)], x='Time (Month-Day/Hour:Minute)', y="Flow from Grid", color= 'Scenario', title="Flow from Grid when system is controlled by the prosumer (in 15-minutes time-steps), rate is TOU", labels={"Flow from Grid":"kWh"}, range_y=[-8,8])
fig.show()
print('The prosumer controls (Scenarios)')
pro_results = pro_results.replace(['s0','s1','s2','s3','s4','s5','Time-of-Use', "Time-of-Use EV & flat", "flat"],['s0: Status Quo', 's1: Battery', 's2: Battery+PV', 's3: EV', 's4: EV+PV', 's5: EV+PV+Battery', "TOU", "TOU EV & Flat", "Flat"])
tmp = pro_results[pro_results.Rate == "TOU"]
Scen = ['s0: Status Quo', 's1: Battery', 's2: Battery+PV', 's3: EV', 's4: EV+PV', 's5: EV+PV+Battery']
for i in range(0,6):
db = tmp[tmp.Scenario==Scen[i]]
fig = go.Figure(data=[go.Sankey(
node = dict(
pad = 15,
thickness = 20,
line = dict(color = "black", width = 0.5),
label = ["Grid (G)", "Solar PV (R)", "ESS (B)", "Load (L)", "EV (E)"]
),
link = dict(
source = [0, 0, 1, 1, 1, 2, 2, 0, 1, 4, 4],
target = [2, 3, 0, 2, 3, 0, 3, 4, 4, 0, 3],
value = [db.loc[:, "x^GB"].sum(), db.loc[:, "x^GL"].sum(), db.loc[:, "x^RG"].sum(), db.loc[:, "x^RB"].sum(), db.loc[:, "x^RL"].sum(), db.loc[:, "x^BG"].sum(), db.loc[:, "x^BL"].sum(), db.loc[:, "x^GE"].sum(), db.loc[:, "x^RE"].sum(), db.loc[:, "x^EG"].sum(), db.loc[:, "x^EL"].sum()]
))])
fig.update_layout(title_text=f"Scenario {Scen[i]}", font_size=10)
fig.show()
pro_key_results = pro_key_results.replace(['s0','s1','s2','s3','s4','s5','Time-of-Use', "Time-of-Use EV & flat", "flat"],['s0: Status Quo', 's1: Battery', 's2: Battery+PV', 's3: EV', 's4: EV+PV', 's5: EV+PV+Battery', "TOU", "TOU EV & Flat", "Flat"])
fig = px.bar(pro_key_results[pro_key_results.Sell == True], x='Rate', y="npv (%6.19)", color = 'Scenario', barmode='group', title="NPV of scenarios where prosumer controls",range_y=[-35000,15000])
fig.show()
fig = px.bar(pro_key_results[pro_key_results.Sell == True], x='Rate', y="GMP's cost", color = 'Scenario', barmode='group', title="GMP's cost when system is controlled by PRO",range_y=[0,6000],labels={"GMP's cost":"System cost ($)","Sell":"Selling rate"})
fig.show()
pro_results = pro_results.replace(['s0','s1','s2','s3','s4','s5','Time-of-Use', "Time-of-Use EV & flat", "flat"],['s0: Status Quo', 's1: Battery', 's2: Battery+PV', 's3: EV', 's4: EV+PV', 's5: EV+PV+Battery', "TOU", "TOU EV & Flat", "Flat"])
tmp = pro_results[pro_results.Sell==True]
fig = px.line(tmp[tmp.Rate == 'TOU'], x='Month', y="Peak load", color = "Scenario", range_y = [0,30], title="Peak loads when system is controlled by (Rate: TOU)", labels={"Peak load":"kW"})
fig.show()
fig = px.line(tmp[tmp.Rate == 'TOU EV & Flat'], x='Month', y="Peak load", color = "Scenario", range_y = [0,30], title="Peak loads when system is controlled by (Rate: TOU EV & Flat)", labels={"Peak load":"kW"})
fig.show()
print('The prosumer controls (Rates)')
tmp = pro_results[pro_results.Scenario == "s3: EV"]
Rates = ['Flat', 'TOU', 'TOU EV & Flat', 'RTLMP', 'TOU&RTLMP']
for i in range(0,5):
db = tmp[tmp.Rate==Rates[i]]
fig = go.Figure(data=[go.Sankey(
node = dict(
pad = 15,
thickness = 20,
line = dict(color = "black", width = 0.5),
label = ["Grid (G)", "Solar PV (R)", "ESS (B)", "Load (L)", "EV (E)"]
),
link = dict(
source = [0, 0, 1, 1, 1, 2, 2, 0, 1, 4, 4],
target = [2, 3, 0, 2, 3, 0, 3, 4, 4, 0, 3],
value = [db.loc[:, "x^GB"].sum(), db.loc[:, "x^GL"].sum(), db.loc[:, "x^RG"].sum(), db.loc[:, "x^RB"].sum(), db.loc[:, "x^RL"].sum(), db.loc[:, "x^BG"].sum(), db.loc[:, "x^BL"].sum(), db.loc[:, "x^GE"].sum(), db.loc[:, "x^RE"].sum(), db.loc[:, "x^EG"].sum(), db.loc[:, "x^EL"].sum()]
))])
fig.update_layout(title_text=f"Rate: {Rates[i]}", font_size=10)
fig.show()
gmpp = gmp_key_results.replace(['s0','s1','s2','s3','s4','s5','Time-of-Use', "Time-of-Use EV & flat", "flat"],['s0: Status Quo', 's1: Battery', 's2: Battery+PV', 's3: EV', 's4: EV+PV', 's5: EV+PV+Battery', "TOU", "TOU EV & Flat", "Flat"])
gmppp = gmpp.drop(["Pro's cost (buy : flat & sell : flat)", "Pro's cost (buy : Time-of-Use EV & sell : Time-of-Use)", "Pro's cost (buy : Time-of-Use EV & sell : flat)", "Pro's cost (buy : Time-of-Use & sell : Time-of-Use)", "Pro's cost (buy : flat & sell : Time-of-Use)", "Pro's cost (buy : Time-of-Use & sell : flat)"], axis=1)
gmpppp = gmppp.drop(["B_max","B_min","U_c","U_d","E_max","E_min","use_pv","price","Agent","Trade","Rate"], axis=1)
gmpppp['Net trade'] = gmpppp["GMP's cost"] - gmpppp["System cost"]
gmpppp['Maximum peak load'] = gmp_results.groupby("Scenario").max()["Peak load"]
gmpppp = gmpppp.rename(columns={"GMP's cost": 'Total cost', "internal rate of return (%)" : "IRR", "saving" : "Saving", "npv (%6.19)": "NPV"})
# gmpppp["Total flow from grid"] = gmp_results.groupby("Scenario").sum()["Flow from Grid"]
tmp = gmp_results.groupby("Scenario").sum()
gmpppp["Total flow taken from grid"] = tmp["x^GL"] + tmp["x^GB"] + tmp["x^GE"]
gmpppp["Total flow sent to grid"] = tmp["x^RG"] + tmp["x^BG * eta_d"] + tmp["x^EG * eta_d"]
temp = ['Scenario', 'Net trade','System cost', 'Total cost', 'Saving','NPV','IRR', 'Maximum peak load', "Total flow taken from grid", "Total flow sent to grid"]
gmpppp = gmpppp[temp]
gmpppp
pro1 = pro_key_results.replace(['s0','s1','s2','s3','s4','s5','Time-of-Use', "Time-of-Use EV & flat", "flat"],['s0: Status Quo', 's1: Battery', 's2: Battery+PV', 's3: EV', 's4: EV+PV', 's5: EV+PV+Battery', "TOU", "TOU EV & Flat", "Flat"])
proo = pro1.drop(["B_max","B_min","U_c","U_d","E_max","E_min","use_pv","price","Agent"], axis=1)
tmp = ['Rate','Sell','Scenario',"Pro's cost","GMP's cost",'saving','npv (%6.19)','internal rate of return (%)']
proo = proo[tmp]
proo.reset_index(drop=True, inplace=True)
pro_results["Total flow taken from grid"] = pro_results["x^GL"] + pro_results["x^GE"] + pro_results["x^GB"]
pro_results["Total flow sent to grid"] = pro_results["x^RG"] + pro_results["x^BG * eta_d"] + pro_results["x^EG * eta_d"]
for s in ['s0: Status Quo', 's1: Battery', 's2: Battery+PV', 's3: EV', 's4: EV+PV', 's5: EV+PV+Battery']:
for p in ['Flat', 'TOU', 'TOU EV & Flat', 'RTLMP', 'TOU&RTLMP']:
for t in [True, False]:
proo.loc[(proo.Sell == t) & (proo.Scenario == s) & (proo.Rate == p), "Total flow sent to grid"] = pro_results.groupby(['Rate','Sell','Scenario']).sum().query(f'Rate == "{p}" and Sell == {t} and Scenario == "{s}"')["Total flow sent to grid"].to_numpy()[0]
proo.loc[(proo.Sell == t) & (proo.Scenario == s) & (proo.Rate == p), "Total flow taken from grid"] = pro_results.groupby(['Rate','Sell','Scenario']).sum().query(f'Rate == "{p}" and Sell == {t} and Scenario == "{s}"')["Total flow taken from grid"].to_numpy()[0]
proo.loc[(proo.Sell == t) & (proo.Scenario == s) & (proo.Rate == p), "Maximum peak load"] = pro_results.groupby(['Rate','Sell','Scenario']).max().query(f'Rate == "{p}" and Sell == {t} and Scenario == "{s}"')["Peak load"].to_numpy()[0]
proo = proo[['Rate', 'Sell', 'Scenario', "Pro's cost", "GMP's cost", 'saving', 'npv (%6.19)', 'internal rate of return (%)', 'Maximum peak load', "Total flow taken from grid", "Total flow sent to grid"]]
# proo.reset_index(drop=True, inplace=True)
pr = proo.query("Sell == True").drop(columns="Sell")
pro_2 = proo.query("Sell == False").drop(columns="Sell")
pr["NPV'"] = pro_2["npv (%6.19)"].values
pr["System Cost'"] = pro_2["GMP's cost"].values
pr.rename({"Pro's cost": 'Total Cost', "GMP's cost": 'System Cost', "saving": "Saving", "npv (%6.19)": "NPV", "internal rate of return (%)": "IRR"}, axis=1, inplace=True)
pr.reset_index(drop=True, inplace=True)
pr.set_index('Rate', inplace=True)
pr